A HYBRID FINITE DIFFERENCE FINITE VOLUME APPROACH 
TO SOLVE FIRST-ORDER HYPERBOLIC CONSERVATION LAWS 
WITH SUPERIOR ACCURACY 
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Abstract. A hybrid finite difference— finite volume (FD-FV) approach for discretization in space 
is proposed to solve first-order hyperbolic conservation laws. Unlike any conventional finite difference 
method (FDM) or finite volume method (FVM), this approach uses both cell-averaged values and 
nodal values as degrees of freedom (DOF). Consequently it is inherently conservative like FVM and 
easy to extend to high-order accuracy in space like FDM. The proposed FD-FV approach works 
for arbitrary flux functions, whether convex or non-convex; and it does not require any exact or 
approximate Riemann solver hence it is also computationally economical. Method of lines is adopted 
for time integration in present work; in particular, explicit Runge-Kutta methods are employed. It 
is theoretically proven and numerically confirmed that in general, the proposed FD-FV methods 
possess superior accuracy than conventional FDM or FVM. Linear stability is studied for general 
FD-FV schemes - both space-accurate and time-stable FD-FV schemes of up to fifth-order accuracy 
in both space and time are presented. Numerical examples show that as long as the solutions are 
smooth, the proposed FD-FV methods are more efficient than conventional FVM of the same order, 
at least when explicit time-integrators are applied. 
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1. Introduction. Hyperbolic conservation laws (HCL) [3 [10] have wide appli- 
cations in areas like computational fluid dynamics (CFD), computational aeroacous- 
tics (CAA), and magnetohydrodynamics (MHD). Since analytical solutions to most 
HCLs are difficult to construct, effective numerical methods become indispensable 
tools for both theoretical analysis and applications. Three desired properties of prac- 
tical numerical schemes for such problems are: (1) conservation, which ensures correct 
shock speed; (2) high-order accuracy when the solution is smooth; and (3) stability. 
In addition, versatility is also preferred, which means that the method should work 
for general flux function that meets minimum mathematical requirements, yet still 
provides correct solutions. Most popular numerical techniques, like finite difference 
methods (FDM), finite volume methods (FVM) and finite element methods (FEM), 
have their own strength and weakness in resolving these issues. This work, on the 
other hand, makes an attempt towards constructing versatile numerical schemes such 
that all these goals can be naturally achieved. 

Finite volume methods, arguably the most popular technique in practice for 
HCLs, is naturally conservative. Usually, these methods are of Godunov type [T2], 
which requires an exact or approximate Riemann solver [281 1111 137] . These Riemann 
solvers often highly depend on the specific flux function, and can be expensive to 
compute, especially when the flux function is non-convex|17j . In addition, since these 
methods typically operate on cell-averaged values, high-order accuracy in space is dif- 
ficult to achieve. Most practical FVMs are only second-order accurate, for example, 
the monotone upstream-centered schemes for conservation laws (MUSCL) [39], flux- 
corrected transport (FCT) [3] and the Kolgan's method [23] • This is partially due to 
the fact that cell-averaged value can be treated as a second-order approximation to 
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the nodal value at cell centroid [211 [20] , However, to construct higher than second- 
order FVMs, one must be very careful to treat the cell- averaging nature of discrete 
data properly, which leads to further complication. Some examples are the piecewise 
parabolic method (PPM) [S], the k-exact method [3 , 27 and essentially non-oscillatory 
(ENO) [17] or weighted essentially non-oscillatory (WENO) [25] schemes. 

On the contrary, finite difference method, being numerical methods based on 
nodal values and conventional Taylor series expansions, are easy to be designed to 
high-order accuracy. In addition, since most discrete operators of EDM are linear, 
von Neumann analysis can be applied for preliminary study of the linear stability 
when designing a scheme. However, these methods often lack inherent conservation 
and additional measures must be taken [TSJ |5J |3S] to ensure the correct shock speed. 

Finally, the third type of methods, EEMs, are not originally designed for advection 
dominated problems, like most HCLs of concern. To account for instability caused 
by advection, two approaches are often employed to enhance or modify the original 
Galerkin type methods. The first approach introduces additional terms [THIET] to sta- 
bilize the whole system, without modifying the Galerkin framework. These methods 
are conservative, as long as constant function is included in the weighting functional 
space |2| . The other approach, represented by the discontinuous Galerkin methods [5] , 
requires numerical fluxes like FVMs at element interfaces, thus they depend on the 
particular physical flux function. 

In this work, a hybrid approach is proposed for spatial discretization, such that 
it combines the advantages of both EDM and EVM in the sense that the method is 
inherently conservative like FVMs as well as easy to extend to high-order accuracy 
like EDMs. In addition, the new approach is suitable for general flux functions, 
convex or not, and it requires no numerical flux functions. To this end, both cell- 
averaged values and nodal values at cell faces are counted as degrees of freedom in 
the proposed ED-FV approach. Method of lines is adopted for time integration; in 
particular, explicit Runge-Kutta time stepping methods are employed in this work 
although it can be extended to implicit methods like conventional EDMs and FVMs. 
For discretizing the cell-averaged values, fluxes at cell faces are computed as exact 
flux functions evaluated at nodal values. To discretize the nodal values in space, 
polynomial reconstruction taking into account of both adjacent nodal values and cell- 
averaged values is performed. This reconstruction leads to a linear discrete differential 
operator (DDO) at cell faces, which enables the von Neumann stability analysis to 
the proposed ED-FV schemes. 

The present paper focuses on ID case. It is proven theoretically and verified nu- 
merically that the proposed ED-FV approach generally has superior accuracy property 
compared to conventional finite difference methods providing the same order of poly- 
nomial reconstruction. This property is demonstrated by numerical example to show 
that it also carry to 2D case on Cartesian grids. Thus the proposed methods are po- 
tentially more computational efficient than conventional numerical methods to solve 
HCLs. 

Some related works in literature are briefly discussed in the following. These 
methods, like the present one, make use of DOEs with different mathematical mean- 
ings as dependent variables. The constrained interpolation profile method (CIP) [42] 
and the interpolated differential operator method (IDO) \T, "W have spatial deriva- 
tives and nodal values as DOEs. They employs Hermite interpolation for local data 
reconstruct. However, these methods are more suitable for elliptic equations and 
parabolic equations in which case the solutions are more regular than HCLs. The 
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enhanced versions to handle conservation laws, naming CIP-CSL |4T] (CSL stands for 
conservative semi-Lagrangian) and IDO-CF (TH] (CF stands for conservation form), 
are able to handle certain hyperbolic equations like the Euler equations and incom- 
pressible Navier-Stokes equations in conservation form. These enhanced versions are 
closely related to the present approach, in the sense that they include cell-averaged 
values into the DOFs and discard all the spatial derivatives. However, compared to 
CIP-CSL and IDO-CF, the current approach is more general and it is suitable for all 
flux functions. In addition, to the best knowledge of the author, provable superior 
accuracy in space is never shown before for alike numerical schemes. 

To this effect, the remainder of the paper is organized as follows. Section [2] 
presents the general framework of the proposed approach for one-dimensional equa- 
tions as well as a simple extension to two-dimensional case. The convergence property 
in space is studied in Section [3] This section also shows that in most cases, the FD- 
FV methods have superior accuracy property than conventional FDMs providing the 
same order polynomial reconstruction. Next, Section [4] studies the linear stability 
property of both the semi-discretized and the fully discretized forms of the proposed 
methods; in the latter case explicit Runge-Kutta time-integrators are employed for 
simplicity. Practical FD-FV schemes that are both high-order accurate and linearly 
stable are presented in Section [5j Numerical examples that confirm the theoretical 
results are offered in Section [6j where a number of benchmark problems are solved to 
show both the accuracy and the robustness of the new approach. Finally, Section [7] 
concludes this paper. 

2. A hybrid finite-difference and finite- volume (FD-FV) approach. Here 
and throughout the remainder of the paper, the following notations are used: 

• The normal fonts {x,y,t,u, etc.) designate scalar variables; and the bolded 
fonts {x,w,f, etc.) designate vector-valued variables. 

• Unless otherwise noted, the subscripts x, y mean partial derivatives in space, 
and the subscript t indicates time derivative. 

• '\{x,t) designates the exact value of a variable f evaluated at the location x 
(in Cartesian coordinates) and time t. 

• The overline f indicates the cell-averaged value of a variable f over a control 
volume C; whereas in absense of the overline, nodal value is indicated. 

2.1. General FD-FV formulation for one-dimensional hyperbolic con- 
servation laws. Consider the system of first-order hyperbolic conservation laws (hy- 
perbolic conservation laws or HCL for short throughout the remainder of this paper) 
in one space dimension: 

wt + f{w), = (2.1) 

where ti; e M"* is a real-valued vector and / : M'* M'' is a differentiable function 
of w. Denote the Jacobian matrix of the flux function as J{w) = df /dw and it is 
supposed to be known analytically. 

Consider the fixed computational domain CI = [a, b] and a uniform mesh with N 
control volumes: fl = ufLiCi, Ci = [a + {i — l)h, a + ih], where h = [h— a)/N is the 
cell size. The cell centers are denoted by xf. Xi ~ a+{i — l/2)h, i = 1, 2, • • • ,N; and 
the cell faces by Xi^i/2- 3^^+1/2 = a + ih, i = 0,1, ■ ■ ■ , N . 

The proposed FD-FV methods seek the numerical approximations to both the nodal 
values at cell faces, w^_^_^/2 ~ iJ''{^i+i/2jt"')i a-nd approximations to the cell-averaged 
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values: 

^ j^J^ w{x,ndx (2.2) 

Here \Ci\ — Xi_^^i/2 ~ = h is the size of the control volume Ci. 

Method of lines is adopted in present work, which means Eq. (2.1 ) is first discretized 

in space to obtain an ODE system for the time-dependent variables: 

Wi+i/2{t) ~ w{x,+i/2,t), w^{t) -r^ 1 w{x,t)dx (2.3) 

JCi 

The resulting ODE system is then solved by any preferred ODE solver. In this paper, 
the explicit Runge-Kutta methods are employed. Here and throughout the remainder 
of the paper, f ^ , t"+i/2 ^^"^ tii ti+1/2 designate approximated cell-averaged and nodal 
solutions to the fully discretized system and semi-discretized system, respectively. In 
the latter case, the dependence in t is omitted for convenience. 

To derive the semi-discretized formula for cell-averaged values, the governing 
equation is integrated over each cell Ci to obtain: 



w{x,t)dx^ ^ 1^ [-^ (■J«(a;i+i/2,<)) - / {w{xi_i/2,t))] = 



di yjc, 

which suggests the following semi-discretized scheme for w 
dwi 1 



dt ■ \C,\ [/(^■'+i/2)-/K-i/2)]-0 (2.4) 



Note that Eq. (2.4) is exact in the sense that there is no truncation error associ- 
ated with this semi-discretization formula. On the other hand, to obtain the semi- 
discretized formula for nodal values, Eq. (2.11 is written in non-conservative form at 
a cell face Xi^i^2- 

■Wt{Xi+i/2,t) + J {w{Xi+i/2, t)) Wa:{Xi+i/2,t) = 

The following numerical scheme is proposed here for solving Wi_^_l/2■ 

'^'^''l^''^ + («^) [2':r'^]i+l/2 =0' = -^^+1/2 (2.5) 

iu does not have to be chosen as ifj_|.]^/2i but it is the simplest and the most natural 
choice and is adopted in this paper. This choice is also shown to simplify a number 
of issues when nonlinear problem is concerned (see Section [5|. is the discrete 
differential operator that characterizes the FD-FV scheme, which takes into account 
both nearby nodal values and cell-averaged values. 

In this paper, linear 2?^ is considered so that von Neumann stability analysis 
applies. It can be expressed in general form as follows: 

\^x'^\+Xl2 = Jl ^I'^^+i + lYl (2.6) 



Here ai and /3/ are constants, and q indicates the stencil of the numerical scheme. 
Note that Eq. (2.6 1 is constructed to a desired order of accuracy based on Taylor 
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series expansions. In particular, supposing sufficiently smooth solution, then for nodal 
values conventional Taylor series expansion leads to: 

°° (Ifi)™- 

w{Xi+i/2+l,t) = w{x,+i/2,t) + ^ —d^w{x,+ i/2, t) (2.7) 

m— 1 

and for cell-averaged values, the following expansion is employed: 

2+lh 



2 /■a:i + i/2+"i 

Wi+i{t) = T I w{x,t)dx 

/ w{x,+i/2,t) + Y\ 'A-^ d^w{x,+T^/2,t) dx 



+ l/2 + {l-l)h 
' [■W{x,+ y2,t) + 

w{x,+,/2,t) + ^ 1 \ [ ^ d:^w{x,+„2.t) (2.8) 

•■^ (m + 1)! 



In theory, Eqs. (2.7 2.8) enable one to construct to arbitrary order of accuracy in 
space. 

Remark: In the special case ai = 0, VZ, Eq. (2.6) is the conventional finite difference 
operators. Thus in this case, the proposed FD-FV methods reduce to classical FDM, 
with cell-averaged values being numerical artifacts. 

2.2. A simple extension to two dimensions on Cartesian meshes. The 

methodology can be easily extended to solve 2D hyperbolic conservation laws in the 
case of Cartesian meshes. The governing equation is: 

wt + f{w),^g{w)y^{) (2.9) 

Let J = df/dw and K = dg/dw be the Jacobian matrices for the fiux functions in 
X— and y— directions, and they are supposed to be known analytically. 
Consider a Cartesian mesh with cells C^.j — [x.1^1/2, a;i+i/2] x [yj_i/2j ?/j"+i/2]j where 
the cell center of C^j is {xi,yj) = {{xi_i/2 +Xi+i/2)/2, {yj-1/2 + yj+1/2)/'^) ■ In the 
semi-discretized form, the cell-averaged dependent variables are naturally defined as: 



1 



'^»j(*)~T7^/ w{x,t)dx (2.10) 



To reuse the DDO proposed in the ID case, the nodal values at the following locations 
are chosen in this work: 

Wij+i/2{t) « w{xi,yj+i/2,t), tOi+i/2j(t) w w{xi+i/2,yj,t) (2.11) 

The locations of the dependent variables are illustrated in Figure |2.1| 



= (2.12) 



The cell-averaged variables are advanced in time as follows: 

dw,^j f {■W,+ i/2,j) - f {Wi-i/2,j) g {w^^j + i/2) ~ g {Wij_if2) 



dt hx hy 



Here and hy are edge lengthes in x- and y-directions of the control volumes. Unlike 
the ID case, Eq. (2.12) is not exact, because edge-averaged fluxes are replaced by flux 
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Fig. 2.1. Locations of dependent variables of the 2D FD-FV scheme on Cartesian mesh: o 
locations of cell-averaged variables; • - locations of nodal variables 



functions evaluated at edge centers, which introduces second-order error in space. 
Higher-order fluxes can be obtained, for example, using nearby nodal values on the 



same axis. However, Eq. (2.12) is adopted for numerical examples in this paper, for 
the reason that it is sufficient to show the superior accuracy property (to be presented 
in Section[3| of the proposed FD-FV schemes in 2D case. FD-FV schemes with higher- 
order fluxes in multiple dimensions will be addressed in future work. 
The semi-discretized formula for nodal values are given by the following: 



2j 



dt 



dw. 



dt 



[i>xW]i j + i/2 + K{w,,j + i/2) [DyW] 



iJ+1/2 







(2.13) 
(2.14) 



Here and Vy are the ID discrete differential operators developed in previous 
section. On the other hand, T)^ and "Dy are classical discrete operators in FDMs that 
only uses nodal values. For example, using first-order one-sided operators, one has 



Forward^ 



1i+l/2j 



>i+l/2,j + l 



or [V. 



Backward^ 



iT^r orward^ .1 



or [V^ 



Backward^ 



']ij + l/2 



''i+lJ+1/2 



'i,j + l/2 



'h,j+l/2j , 



'h^l,] + l/2} 



Finally, before going into analysis of the accuracy and stability properties of the pro- 
posed method, a few remarks should be addressed: 

Remark 1: This section provides the general framework of the FD-FV schemes. To 
construct stable and accurate FD-FV methods to solve hyperbolic problems, adequate 
upwinding property (or more precisely, linear stability) must be considered. These 
are provided in Section [5] where the practical FD-FV schemes are presented. 
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Remark 2: The 2D FD-FV scheme framework presented is a simple one. Neverthe- 
less, it is sufficient to illustrate the superior accuracy property (see Section [s]) in 2D 
case by this simple framework as shown in Section [6] 

Remark 3: In the 2D case, locations of dependent variables other than the one in 



Figure 2.1 can be explored. For example, one can choose (a;i+i/2, yj+1/2) to be the 
points where nodal values are evaluated. A discussion of the advantages and disad- 
vantages of each of these configurations is outside the scope of this paper, and will be 
addressed in future work. 

Remark 4: Some related works |1TJ[T2] introduce edge-averaged values along control 
volume edges as dependent variables. However, since the flux function and edge- 
averaging operator typically don't commute, one still needs to: (1) approximate edge- 
averaged fluxes with edge-averaged conservative variables; and (2) develop governing 
equations for the edge-averaged variables. To avoid such complications, edge-averaged 
variables are not considered in this work. 



3. Accuracy analysis. Consider applying the scheme given by Eqs. (2.4 2.51 
to solve ID scalar advection equation with constant propagation velocity: 

Ut -I- cux — (3.1) 

The flux function is given by f{u) — cu, and the Jacobian is J(u) = c. The accu- 
racy and stability properties of the proposed FD-FV scheme are studied by applying 



Eqs. (2.4 2.5) to solve Eq. (3.1) with simple wave solution 



To this effect, the discrete operator that characterizes the FD-FV scheme. 



which is given by Eq. (2.6 1, is divided into two operators A and Af as follows: 

^ q ^ q 

l=-q+l l=-q 



And it follows: 



-1/2 - L"^"^Jj+l/2 ^ l'^ '^Ji+1/2 

Note that no subscript x is introduced to A or Af, because they do not necessarily 
constitute a consistent discrete differential operator that approximates first-order spa- 
tial derivative at And it is shown later that, this is the key to the superior 
accuracy property of the proposed FD-FV schemes. 

3.1. Application of the proposed schemes to ID advection equation 



with simple wave solutions. Consider the simple wave solution to Eq. (3.11 with 
wave number k: 

u(a;,t) = e^'=(^-'='' (3.3) 



In the following, the superscript ★ designates the exact values corresponding to (3.3). 
To avoid confusion, i is reserved for imaginary unit in the remainder of this section, 
and j will be used as a generic subscript that indicates a particular control volume. 
It is also supposed that the mesh is uniform with cell size h, and Xj — jh, Xj^i/2 = 
(j + l/2)h. Then M*(t) and w*_^_-^y2(0: ^^'^ exact cell-averaged value and nodal value 
at time t, are given by: 



1 



U + l/2)h 



^H^-ct)^^ = le--'=*e^J« (e''/^ - e-''/^ (3.4) 



'^3= T. ^ 
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Here 9 — kh. Define the following quantities for convenience: 

A*{t) = N*{t) = e-"'''\ RAO) = -e''<^ ( e'^/^ - e"^/A , 

id \ / 

Then the exact solutions (3.4-3.5) are: 



S. 



J + l/2 



{0) 



„j(j+i/2)e 



(3.6) 



Next apply the FD-FV scheme ( plp^ to solve Eq. It follows that the ODE 

systems obtained from discretization in space for dependent variables Uj and Uj+1/2 
are: 



du 



du 



'3 + 1/2 

dt 



(3.7) 
(3.8) 



Here A and N are given by Eqs. (3.2 1. 



To find the solutions to the ODE systems (3.7- 3.8[ ) with exact initial data: 

It is supposed first and verified afterwards that the solutions are given by: 

= A{t)R,{e), 7.,+i/2 = N{t)S^+y2{0) 



(3.9) 



with A{0) = iV(0) = 1. Note that both A{t) and N{t) also depend on 9. However, 
to simplify the notations and also to be consistent with the vocabulary "ODE" that 
follows, the dependence on 9 is omitted but should be kept in mind. Then it follows 
that A{-) and N{-), which are approximations to A*{-) and N*{-), satisfy the ODEs: 



A\t)+ickN{t) = 
N'{t) + cka{9)A{t) + ckb{9)N{t) = 



Here a(-) and b{-) are functions of 9: 



1 - e- 



(3.10) 
(3.11) 



(3.12) 



l=-q+l 



l = -~q 



The ODEs ( |3.10||3.11| are obtained by inserting Eqs. (|3^ to Eqs. ( |3.7p.8[ ) and 
applying the following equalities: 



i9Rj{9) — Sj+i/2{9) - Sj_i/2{9), Sj+i/2+i{9) = Sj+i/2{9)t 
RMe)/S,+i/2{0)^^,e^'\l-e-'') 



iW 



lytical solution is given by: 



The ODE system (3.10-3.11) is the key for accuracy and stability analysis. Its ana- 



A{t) 
N{t) 



_ ^-ckc{e)t 
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The eigenvalue decomposition of the coefficient matrix C{9) is: 

A2 1 



C{9) 



i i 




Ai 0' 








A2 





i(A2-Ai) A2-A1 
Ai 1 
i(A2-Ai) A2-A1 



(3.14) 



where the eigenvalues are functions of 9: 

Ai(^) - I (b{9) + ^b{9)^ + 4ia{9)') , A2W - ^ [b{9) - ^ b{9)^ + Ata{9)'^ (3.15) 
Consequently, the solutions are expressed analytically as: 



Ait) 



A2 — Ai 



-ckXit 



-ckX^t 



A2 — Ai 



N(i\ - -^l(-^2 -^) cfcAit _ M^l^lilp-cfcA2t 
*(A2-Ai) *(A2-Ai)' 



(3.16) 
(3.17) 



Note that the dependence on 9 in Ai^2 is omitted here and throughout the remainder 
of this section for convenience. 



Since Rj{9) and Sj^i/2{9) are known, it is easy to verify that Eqs. (3.16 
to the exact solutions to the semi-discretized problem (3.7 3.8) through 



3.17) lead 
^:qs:^(|3^. 



Comparing with Eqs. (3.4-3.5), the order of spatial accuracy of the FD-FV scheme, 
or the expected convergence rate in space, is defined as follows: 

Definition 3.1. The FD-FV scheme defined by the DDOVx ispth-order (p>l) 
accurate in space, if: 

A^^^l + Cait)9P + 0{9^+'), ^^=l + c„{t)9P + Oi9P+') (3.18) 

Here Ca and c„ are functions of time that are independent of 9, and at least one of 
them is not zero. 



Comparing the exponential terms in Eqs. (3.16-3.17) and noticing that A*{t) = 
N*{t) — e^*'^'"'*, one of the two eigenvalues Ai,2 is a numerical approximation to 
the imaginary unit i, whereas the other one is spurious. Without loss of generality, 
Ai is chosen to be close to i and A2 is supposed to be spurious. 
The designed order of accuracy of the DDO is defined as follows: 



Definition 3.2. The DDO given by Eq. (2.6) is pth- order (p > 1) accurate 
in space, if: 



['D..v%+i/2 = d.v*+y2 + CpdP+'v*^,/2h^ + 0{hP+^) 



(3.19) 



for any sufficiently smooth function v*{x,t). Here Cp is a non-zero constant that is 
independent of the function v* and cell size h. 

In the case of conventional EDMs (or A = as explained before), the numerical 
scheme usually has the same order of accuracy as V^. However, it is to be shown next 
that this is not the case for most FD-EV methods proposed in this paper. 
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3.2. Accuracy analysis of and the corresponding FD-FV scheme. In 

order to study the accuracy of the FD-FV scheme that is determined by the operator 
Vx, the relationship between the coefficients a;, /?; and the order of is first estab- 
lished by the following proposition: 



Proposition 3.3. The DDO given by Eq. {2.6) is pth-order (p> I) accurate 
with leading error coefficient Cp ^ if and only if the constants ai, —q + 1 < I < q 
and Pi, —q <l<q satisfy: 

a; + ^ A = (3.20a) 

E -^ai+Y.iPi = l (3.20b) 

E ( \J A=0, m = 2,-..,p (3.20c) 

^-^ (m + 1)! ^-^ m! 

l = -q+l ^ ' l = -q 

1 lp+2 _ (J _ T\p+2 1 jp+1 



The proof is straightforward: the Taylor series expansions (2.7 2.8) of cell-averaged 
variables and nodal variables are plugged into Eq. (2.6), then Eqs. (3.20) are obtained 



according to Eq. (3.19). 
Defining the following quantities 
^ q 

Cl'rn. 



(to + 1) 



1 (r'+i-(/-l)'"+i)a,, 6,„ = J-VrA, TO = 0,1, 

-f 1)! to! 



l=~q+l 



i=-q 



then Eqs. (3.201 are simplified to: 

ao + &o = 0, ai+6i = l, a,„ + = 0, to = 2, • • • ,p ap+i+hp+i=Cp (3.21) 

In addition, from Eqs. (3.12), one obtains by Taylor series expansions of e*'^ the 
following expressions of a{9) and b{9) in terms of and bm'- 



m=0 

oo 



;m + ll 



(3.22) 



(3.23) 



Suppose is designed to be pth-order accurate in space, so that Eqs. (3.21 ) hold, it 
follows: 

oo 

a{e) + b{e) = I + r(0)0P, r{e) = *"+^+'(am+p+i + 6™+p+i)e"^ 

T?l = 



To compute the two eigenvalues of C{9) using Eqs. (3.15), the term in the square root 
operation is first evaluated: 

2 



¥ + 4m = b^ + Ai{i -b + rOP) = {2i - by + AirO^ =\2i-b 



2ireP 
2i-b 



{2i - 6)2 
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Note that it is assumed here that b{9) ^ 2i, which happens in the special case 60 = 
0, 61 = 2, bm = 0, m > 2 according to Eq. (3.23 1. The last term in the parentheses 



has a crucial impact on the accuracy of the FD-FV scheme, and it is expressed as 



follows by applying Eq. (3.23) and the leading terms of r{6): 



2irep _ 2iP+^cpep+^ + o{ep+^) 



2i -b 2i9- YZ^Q ^"'b. 
To study the FD-FV scheme's spatial order of accuracy, the index is defined: 

Definition 3.4. Given the operator J\f, its index sp is defined by: 




(3.24) 



bo^O 

sp= { I &o = 0, &i ^ 2 

s>2 &o = 0, 61 = 2, 62 = • • • = bs-i = 0, bs^O 



(3.25) 



Consequently, the leading term of the denominator of Eq. (3.241 is of order 0{9^^). 
Hence one establishes: 



2ir6 



2i-b 2x{sf,=i} - b. 



(3.26) 



S/3 



Here Xfs^i} takes the value 1 if = 1 and otherwise. Thus the square root in 
Eqs. (ISjII) is: 



v/&2 + 4za = 2i-b+ E_0P+i"^P + 0(6lP+2-'',8 ) + 0(61^) (3.27) 

2x{sf<=i} - bs„ 

The last term accounts for the contribution of fZ ^i^i ~ (9(6l2p+2-2s^^ ^^le expansion 



of b^ + 4ia. In particular, suppose 2i — b 



(2i-fc)2 



gp+i-sf< ^ O{0^), one has: 



mm{p + q,2p)^2p + 2~ 23/3 
Note that the particular choice of the square root from the two available values in 



Eq. (3.271 is due to the consideration that Ai should approximate z, see the discus- 



Eqs. (3.15) 



sion following Definition 3.1 To this end, by applying Eq. (3.27), one obtains from 



A2 = ^ - i 



2x{s3=i} - bsf, 

00 



(3.28) 



2x. 



iP_5/P+l-S3 ^ Qj-^min(p+2-s^,p)-j (-3 29) 



Note that Ai at best approximates i to (p + 1 — s^)th-order accuracy in space. Thus 
to maximize the accuracy of Ai, should be zero, that is, 60 7^ 0. This leads to the 
major result in this paper: 



Theorem 3.5. If sp = 0, or equivalently bo — /?/ ^ 0, a pth-order (p > 1) 

l=-q 

accurate operator leads to a (p + l)th-order accurate FD-FV scheme provided also 
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that cbo > 0. 



Note that the Definitions |3.1H3.2| are used here. They do not guarantee convergence 
of the underlying FD-FV schemes unless stability is also addressed, which is provided 
later. 



Proof. Given sp — 0, there is X{sf,=i} — and consequently: 



2i-b+ ^0P+i-^/3 

2X{s3=i} - bs„ 

\ m=0 / ° 



= 0{9-^) 



Thus (J = — 1, and p satisfies min(p— 1, 2p) = 2p + 2, or equivalently p = 2p + 3. Then 



Eqs. (3.28 -3.291 reduce to: 



Xi=i-'-^9P+' +OieP+^), A2 = ^ + (6l-l)^ + O(0) (3.30) 



Noticing that: 



A2 - Ai = + Aia = -2i + b + 



and following Eqs. (3.16 3.171 one has: 



A*{t) 



A2 — Ai 



A2 — Ai 



N*{t) { ^ i{\2-\i)) ^{X2-Xl) 



(3.31) 
(3.32) 



The leading terms of every single entity are evaluated: 
Xi-i _ -'{iP+^Cp/bn)0P+^ + 0{eP+^) _ iP+^Cp 



A2 — Ai 
A2(Ai-i) fbo 



bo/0 + 0(1) 



/)2 
'^0 



9P+^ + O{0P+^) 



«(A2 - Ai) \i 



0(1) 



np+2 



o{ep+^) 



p+i 



bo 



+ 0{9P+^) 



bo 

^-ckiX2-i)t ^ ^-ckt[bo/e+(bl-2)z+0(e)] ^ ^-Ehll^-rck(bi~2)t ^ q^q^-^ 



Since h 

e 



> 



0, one has e^'^*'"*/'* < 1 due to the assumption that cbo > 0. Hence 
ck{\2-i)t _ o[l) and it goes to zero exponentially as /i — )• 0. Thus the spurious 
modes associated with A2 can be safely removed from Eqs. ( |3.31 3.32 ) for polynomial 
order-of-accuracy analysis, and they are expressed as e{9) in the following formula. 
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Finally, it is obtained: 



N*[t) \ bo J \ bo J 



bo 



which completes the proof that the corresponding FD-FV scheme is {p + l)th-order 
accurate in space according to Definition |3.1| □ 



Remark 1: The key to the superior accuracy property demonstrated by this 

q 

theorem is — or fog = /3; ^ 0. This explains why it can not be expected for 

l=-q 

conventional finite difference schemes: FDMs correspond to that satisfies A = Q, 
q 

in which case oq = ^ a; = 0. Thus the consistency condition (the first equation 

l=-q+l 



in Eqs. (3.21 )) enforces bo = 0, hence sp > and the order-of-accuracy of the scheme 
is at best the designed order-of-accuracy of the discrete differential operator. 

Remark 2: In the case Sjs = 0, the spurious modes associated with A2 has a 
factor e~'=''«*/'\ which can either decay or grow exponentially as the mesh is refined, 
depending on the sign of cbo- Thus it is crucial to design such that bo has the 
same sign as c, the propagation velocity. This suggests that for conservation law 
systems, the characteristic form of the governing equation should be employed, which 
requires the system to be hyperbolic. Note that the condition cbo > can be regarded 
as a mathematical extension (or explanation) of the "upwind" concept for solving 
advection dominated problems. 

Remark 3: The motivation to use FD-FV schemes that has = is not solely 
due to the proven superior accuracy property in space. A potential danger of designing 



FD-FV scheme with such that oq = 60 = is: in this case, Eqs. (2.4 2.51 allow 
the "constant" solutions in the form Wj = Wa, Vj, Wj_^_i^2 — '^b, where Wa and 
Wi, are arbitrary constants. Since Oq = &o = 0, there is no mechanism in the numerical 
scheme that enforces Wa = w^. Thus in the case sp 7^ 0, the FD-FV scheme allows 
unphysical constant solutions. 

Throughout the remainder of this paper, it is always assumed that 60 7^ 0; hence 



Theorem 3.5 applies. Several such applied at a;j+i/2 with up to fourth-order 
accuracy m space are listed in Table |3.1| In the table, "F-biased" and "B-biased" 
stands for forward-biased and backward-biased respectively. 

4. Von Neumann Stability analysis. This section studies the linear stability 
of the proposed FD-FV schemes in the case bo ^ 0. To this end, the stability of the 



solutions to the ODE systems (3.7 3.8) due to discretization in space is first studied 
in Section [41] Next, time-integrators are considered in Section [T2I where asymptotic 
behavior when /i — is studied. This asymptotic analysis only requires the knowledge 
of 60 of the spatial operator Vx- Finally, sample fully discretized FD-FV schemes are 
studied in Section l473l 
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Order 


[25a:ti]j + l/2 


1st 


Forward: 2 
Backward: 2 


^Mj + l - Uj + 1/2) /ft 
Mj + l/2-Uj)/ft 


2nd 


Forward: 2 
Backward: 2 


'-Mj+3/2 + 3u^- + i - 2Uj + i/2) /ft 

,2uj+i/2 - 3wj + Wj-1/2) /ft 


3rd 


Forward: 
F-biased: 
B-biased: 
Backward: 


- + 17Uj+i - lOUj + 1/2) /(2/l) 
-2Uj+3/2 + lUj+l - 4Uj + i/2 - Uj) /(2ft.) 

lij+i + 4uj+i/2 - 7ltj + 2uj_i/2) /(2ft) 
lOuj+1/2 - lluj + 8uj_3/2 - Uj-i) /(2ft) 


4th 


Forward: 
F-biased: 
B-biased: 
Backward: 


^-2mj-+5/2 7uj+2 - 16uj+3/2 + 23uj+i - 12uj+i/2) /(2ft) 
'mj+2 - 12uj+3/2 -1- 31uj+i - 18uj+i/2 - 2uj| /(6ft) 
2?Ij+i + 18uj+i/2 — 31uj + 12uj_i/2 — /(6ft) 
^12uj+i/2 - 23uj -1- 16uj_i/2 - 7uj_i + 2uj_3/2) /(2ft) 



Table 3.1 
Sample "Dx with bo ^ 



4.1. Stability analysis of the semi-discretized form in the case 60 7^ 0- 

Suppose c > 0, then the "forward" and "backward" directions in Table 3.1 correspond 

to the "downwind" and "upwind" directions respectively. 

Looking at Eqs. (3.16 3.17), there are five sources of numerical errors: 

• Dispersion: the computed wave number is fc = Re{kXi/i); or more conve- 
niently one may study 6 = Re{Xi9/i). 6 is related to 9 by the first equation 
in Eqs. (3.30). Their difference represents the error in propagation speed of 
the simple wave. 

• Dissipation: the exact solution does not damp or amplify in time; but the 
numerical solution may do. This error is measured by £ = Im{Xi9 /i). Given 
c > 0, a negative e is desired since in this case the magnitude of the numerical 
solution does not blow up as time increases. 

• Stationary dispersion: this phase error is associated with t he argume nts of 
time-independent quantities in front of g-<^K^i-^)'t Eqs. (3.16 3.17). The 



corresponding errors in the cell-averaged values and nodal values are: 



= arg 



A2 - i 
A2 — Ai 



<y3a = arg 



Al(A2-^) 
z(A2-Ai) 



(4.1) 



Note that the effect of these terms in the solution is independent of time, but 
only depends on the mesh size. 

Stationary dissipation: this error is associated with the magnitude of the 
same quantities above. And they are defined by: 



Vn 



Xo - i 



X2 — Ai 



Ai(A2-i) 



»(A2-Ai) 



(4.2) 



Again, these quantities only affect the solution via mesh sizes, no matter how 
long the ODEs are integrated. 

Noise: unlike most EDMs, the proposed FD-FV methods have a spurious 
mode in the solution and it is called "noise" in this paper. This quantity is 
measured by the magnitude of the second terms in Eqs. (3.16 3.17). 



^1 ~ ^ C-cfc(A2- 

A2 — Ai 



A2(Al-^) _ 
*(A2-Ai)'' 



ck{X2—i)t 



(4.3) 
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Fig. 4.1. Dispersion (6, left) and dissipation (e, right) properties of selected T>x 



These quantities are more difficult to evaluate since they are functions of both 
wave number and time. However, as described before, A2 is a function of 
space only and the leading term is bo/O; thus providing cbo > 0, as required 
by Theorem |3.5[ the exponential terms in V and v are expected to decay 
exponentially with increasing t and deceasing 9. In the plots provided later 
for selected in Table 3.1 the time is fixed to t = h/c, which corresponds 
to one time step size computed with Courant number 1.0. 



Next, these errors are plotted against 6 for the operators listed in Table [3J^ Note, how- 
ever, only the with Bq > are chosen due to Theorem |3.5| and the assumption that 
c > 0. In particular, the chosen operators are the following: (Ist-order back- 

ward), P^.^P (2nd-order backward), 2?3,«p-&»ased (3rd-order B-biased), Pi^-^P (3rd- 
order backward), l?4,«p-fc»ased (4th-order B-biased) and 2?;^^"^ (4th-order backward). 

the 



4.1 



The dispersion relations between 6 and 6 are provided in the left of Figures 

The four stationary quantities (p^^, 
The noises V and v at 



dissipation is plotted in the right of Figures 4.1 
^„ and ifn are plotted as functions of 



(fa 



in Figures 4.2 



fixed time t = h/c are plotted as functions of 9 in Figures 4.3 



Note that the 0-abscissa is chosen to be in the range [0, tt], in which case the largest 
value corresponds to one cell (or two DOFs) per wavelength. The left figure in Fig- 
ures |4T] illustrates: high-order is more accurate than low-order ones; and given the 
order p, the upwind-biased version is more accurate than the upwind version. How- 
ever, looking at the e-9 plot in the right part of Figures |4.1[ it is clear that the fully 
upwind 4th-order accurate V^^'^p is unstable since e > when 9 tt; whereas the 
4th-order upwind-biased operator 2?4,jip-bia5ed j-^^^g such stability problem. This 
observation is against common intuition that "upwind" is closely related to positive 
numerical viscosity and hence it tends to stabilize the numerical scheme. Thus, one 
must pay special attention not to equalize "upwind" property with linear stability 
when designing FD-FV schemes, especially when high-order accuracy in space is con- 
cerned. Figures |4.2| indicate that the stationary errors do not blow up for the range 
of 9 of interests and for all Vj; considered. Figures [473] show that the noises due to the 
spurious mode are bounded and close to zero in the range of 9 for all ; whereas it 
is again observed that V^-'^p, the 4th-order fully upwind operator leads to the largest 
noise among all six operators. 

4.2. Asymptotic stability analysis of fully-discretized problem as h ^ 0. 



The ODE system (3.7 3.8) must be supplemented with a time-integrator in order to 
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Fig. 4.2. Stationary quantities: (1) upper-left, dispersion for cell-averaged values (ifia)> ('^) 
upper-right, dispersion for nodal values (fajj (3) lower-left, dissipation for cell-averaged values 
('fin)' (4) lower-right, dissipation for nodal values (fn) 



Exatt 

Oi'"" 

^p;-"p 



Exact 

pi.^p 

_e_T}^ "P 
^.-pj "P 




Fig. 4.3. Noises associated with spurious mode X2: (1) left, for cell-averaged values (v); (2) 
right, for nodal values (v) 



construct a practical numerical solver. In this section, several explicit Runge-Kutta 
methods up to fifth-order accurate in time are chosen to illustrate the linear stability 
property of selected DDO ■ 

The building block of these time-integrators is the forward Euler (FE) method, which 
is provided below for the general formula (2.4 2.5) given fixed time step size /S.t = 



At 



fi^l+m) /K-1/2) 



iit"+i — 111" 

%- + l/2 - ^3 + 1/2 



Jj+1/2 
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RK2 


Ki = /C(T^") 






K2 = IC{W" -f- 






VK"+^ = VK" 4 


- iA<(Ki +X2) 


RK3 


Ki = IC{W") 






K2 = /C(VF" + 


AtKi) 




K3 = ICiW" -f- 


lAtK, + \AtK2) 




VK"+^ = VF" 4 




RK4 


Xi = /C(W") 






^2 = /C(W" + 


lAtKi) 




= /C(VF" + 


lAtK2) 




^4 = /C(W" + 


AtKs) 






-\At{Ki+2K2 + 2K:, + Ki) 


RK5 


Ki = /C(VF") 








AtKi) 




^3 = /C(VF" -f- 


\AtKi + \AtK2) 




K4 = /C(VF" -f- 


^AtKi + iiAtK2 - iiAtKs) 




K„ = /C(VF" - 


I'AtKi - lAtK2 + j^AiKg + lAtKi) 




Kg = /C(W" - 


l^AtK2 + iiAtK:, + \AtK^ + l^AtK^) 




VK"+^ = VK" 4 


- ^At {7Ki + 7K3 + 32Ki + 12K5 + 32Kg) 



Table 4.1 



Explicit Runge-Kutta methods with up to fifth-order accuracy 



Simplify the formula by defining the operator /C as follows: 
^n+i ^ + AtKi, Ki ^ ICiW^ 



(4.4) 



Here W is the collection of variables Wj and Wj^i/2 



Then the considered Runge- 
Kutta methods [13l ESI [16] are summarized in Table 4.1 Note that the particular 
Runge-Kutta methods are chosen for the purpose of matching the expected order of 
accuracy of the FD-FV schemes. The schemes in Table |4T] are chosen for simplicity in 
presentation due to the simple coefficients. Optimized schemes, for example, towards 
stability preserving [TH [351 [13 or low storage [101 (HI [53] considerations, are not 
explored here. 

Before combining with Dj; in Table |3.1[ asymptotic stability property of each time- 
integrator in Table [4T] as /i — >■ 0, or equivalently — > 0, is studied. 



To this end, first consider applying FE method to solve Eqs. (3.10-3.11 1, which leads 
to: 

^n+i ^ _ ixON", iV"+i = A^" - A6' {a{e)A" + b{0)N") 
Here A = cAt/h is the signed Courant number. Employing the coefficient matrix C{9) 



given in Eqs. (3.13), the equations above become: 

= {I-X0C{e)) 



]\[n+l 



A" 
TV" 



Providing eigenvalue decomposition of C{9) by Eq. (3.14) and initial condition A" 
N'^ — 1, one has: 







i i 




■(I-A6IA1)" 


iV" 




Ai A2 




(1 - A6IA2)" 



i(A2-Ai) 
Ai 



A2 



i(A2-Ai) A2-A1 





'1 




1 



(4.5) 
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FE 


RK2 


RK3 


RK4 


RK5 


|A|< 


2/I&0I 


2/I60I 


2.51/ |5o| 


2.78/ |6o| 


2.52/ |6o| 



Table 4.2 

Necessary condition for linear stability obtained by asymptotic analysis as h ^ 



Order of FD-FV scheme 


Time- integrator 




bo 




^max 


2nd 


RK2 




2 


1.0 


1.0 


3rd 


RK3 


^x 


6 


0.418 


0.409 


4th 


RK4 


, up— biased 

^ X 


3 


0.926 


0.808 


4th 


RK4 


^X 


9 


0.309 


0.309 


5th 


RK5 


-j^A. up— biased 

^ X 


5 


0.504 


0.494 



Table 4.3 



Full FD-FV schemes combining "D^ with proper time-integrators 



Thus a necessary condition for the method to be stable is: 



11- A^Ai,! < 1 



(4.6) 



Recalling Eqs. (3.30 1, in the limit 6* — >■ 0, a necessary condition of the Courant number 



for achieving linearly stability is given by: 



lAI < 



2 

N 



Similar analysis can be conducted for other methods in Table [4?l] and the correspond 



ing bounds are reported in Table 4.2 



Though not mathematically proven, it turns out that for the combinations of 'Dx and 
time-integrators considered in this paper, Table [42] provides convenient estimates of 
the largest allowable Courant number for the given FD-FV scheme to be stable. 

4.3. Stability property of various FD-FV schemes. The spatial operators 



studied in Section 4.1 are combined with proper time-integrators provided in Sec- 
tion |4.2| to give full FD-FV scheme s. T he chosen combinations are provided in the 
second and third columns in Table 4.3 Note that 2?^, up jg j-^q^ included, since it is 
shown to be unstable on the semi-discretized level in Section |4?T] Also note that the 



time-integrator is chosen to be one-order higher than the designed order of the pairing 



I? 2:- This is due to the superior accuracy property proved in Theorem 3.5 



The solution to the fully discretized scheme providing simple wave initial condition 
can be expressed as: 



A" 

TV" 



p{\ec(o)) 



see for example Eq. (4.5) in the case of forward Euler time-integrator. Here P(-) is 
a polynomial that is solely determined by the time-integrator. The FD-FV scheme is 
stable if and only if: 

|P(A0Ai,2(0))| < 1 (4.7) 

Thus to illustrate the stability property of each FD-FV methods listed in Table |4.3[ 
the contours of maximum value max (|P (A6'Ai(6'))| , \P {\9\2{0))\) is plotted on the 
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Fig. 4.4. Contours o/ max (P(AeA,(e))) on 9 - \ plane: (1) upper-left, V^''"'' -RK2; (2) upper- 
right, Vl'^'^-RKS; (3) middle-left, vl''^^'^^"'""^ -RK4; (4) middle-left, Vl'^P-RK4; (5) bottom, 



9 — X plane, see Figures |4.4| In the figures, as 9 0, the contours with value 1 
intersect the A- axis at the values Ag-s-o, which are provided in Section [4.2| by asymp- 
totic analysis. These values are reported in the fifth column in Table |4.3l The last 
column provides the maximum allowable Courant number Xmax for linear stability 
of the FD-FV scheme. It is clear by comparing the last two columns that, Xe^o is 
a good estimate to X^ax, a-nd the computation of the former value is much simpler 
than the latter. 



5. Accurate and stable FD-FV schemes in practice. The FD-FV scheme 
provided by Eqs. ( 2.4 2.5 1 needs to be polished for practical use. To this end, at least 
two issues must be addressed: 

• First, the stability analysis in previous section indicates that the operator Vx 
must be designed to use data in a certain stencil which takes into account of 
the direction of characteristic velocity. This is not refiected in Eqs. (2.4 2.5). 

• Second, boundary conditions must be treated properly. 
These issues are addressed in following sub-sections. 



5.1. Applying DDO to characteristics. Due to stability consideration, 
should be applied to characteristics rather than components of w, thus Eq. (2.51 
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should be modified, with the aid of the eigenvalue decomposition oi J — J {w): 

J = Rki 



Here A is a real diagonal matrix with diagonal elements Ai < A2 < 



< A<j. R 



[fi 7-2 - ■ ■ Vci] and L = R 
eigenvectors associated with Ai 



where r^jij e M'* are right- and left- 



1, 2, • • • d. Consequently Eq. (2.51 can be modified 



as: 



dw 



i+1/2 



dt 



i+l/2 







(5.1) 



Note that this is a numerical approximation to the linearized equation of Eq. ( |2.1[ ) : 

wt + Jwx = (5.2) 

With the choice w = (see Section[2j), the equation above is exact; in this case, 

J(iOi+i/2) is denoted by Ji+1/2 for short. If w is selected differently for a pth-order 
accurate FD-FV scheme, J must be constructed carefully such that J = Ji+i/2 + 
0{hP); usually, this means that w should be chosen such that w = i«i+i/2 + 0{h'P) 
given sufficiently smooth flux function /(■). This issue does not arise with the choice 
in this paper. 



Eq. (5.1) leads locally to a system of decoupled advection equations: 

duj''.' 



dt 



0, 



(5.3) 



Then the chosen D^. in Section |4] are applied to compute approximations to space 
derivatives of each cj'". For example, if 254, up-biased ci^osen, the scheme described 
by Eq. (5.3 1 is as follows: 

[^2:'^"]i+l/2 = 



\vy- 



i+1/2 



U! 



i+l/2 



Am > 
Ar„ < 



Here the meaning of the notations in the superscripts of are given in Table |3.1 



5.2. Boundary condition treatments. Implementation of Dirichlet type and 
Neumann type boundary conditions are considered here for pth-order accurate FD-FV 
scheme and for forward Euler time-stepping method. For simplicity, one- dimensional 
scalar problem is supposed: 



+ f{'^)x = 



(5.4) 



In addition, the chosen DDO is denoted by T>p~^, which is designed to be {p — l)th- 
order accurate. Without loss of generality, suppose the boundary condition at xi^2 = 
0, which is the left end point of the computational domain, is given by: 



u;(0,i) - .9d(<), 
or uj.j;{0,t) = gn{t), 



Dirichlet EC 
Neumann BC 



(5.5) 
(5.6) 



First consider the Dirichlet boundary condition (5.5). To enforce this boundary condi- 
tion at such that the solution at next time step t"^-^ can be obtained, the following 
steps are taken: 
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1. Set ujIi^ = .9d(<"). 

2. If for some i > 0, [P^^^o;"] ^.^-|^^2 requires inputs w"^-^y27 J < that arc not 
available, a different 2?^ i.'+i/2 ^j^j-^. j^^) a biased stencil that do not require 
any w"_^-^y2 7 3 < 0; (2) has fog 7^ and is designed to be (p — l)th-order 
accurate in space, is applied at Xj_|_i/2 instead of 2?^"^. 

For example, if pi.Sacfe^ard^ p2,Bacfe™ard jj3,B-Mased applied around X = 
in the FD-FV scheme, only the first step is required since none of these operators 
requires information more than one cell away. On the other hand, if D^'^"''^^'^"''^'^ or 
jyi^,B-Uased ghogen instead, the operator at 0:1+1/2 should be modified to 2?3,B-&»ased 
and 2?4.^"-''«ased respectively. 



Next, consider the Neumann type boundary condition (5.6 1. It is enforced at 
as follows: 

1. Replace oj^l^^ with the value such that \pPr-^-Forward^n-]^^^^ ^ ^^^^^n) jg gj^^_ 
isfied. Here 2?p-i,^'o'-'"'fd ig the DDO that is designed to be (p - l)th-order 
accurate, with fog 7^ and only use variables in the forward direction. 

2. If for some i > 0, ^_|_j^^2 requires non-existent variables, the treat- 
ment is the same as in the Dirichlet boundary conditon case. 

The first step can be illustrated by considering, for example, the 4-th order accurate 
FD-FV scheme {p = 4) with operator vl^Backward^ r^^^^ ^ ^-^^ ^^^j^g ig 

computed, such that: 



-V^x ^ \ 1/2 - ffnl^ 



or equivalently 



= ^ (^2 - 8^^3/2 + 17^' - 2/ig„(i")) 



Remark: In practice, to solve first-order equations like the HCL considered in this 
paper, the spatial order of accuracy of the numerical scheme can be one-order lower 
at the boundary without affecting the overall order of accuracy in space [TS]. Thus 
the previous implementations of boundary conditions can be simplified by employing 
one-order lower TD^ at and/or near the boundary. 

6. Numerical examples. The proposed FD-FV schemes in Table [43| are tested 
by solving various ID problems. The superior accuracy property is demonstrated 
by solving ID advection equations and ID Euler equations with smooth solutions. 
Transport problems with discontinuous solutions are also tested, in which case linear 
stability analysis does not apply. As expect, Gibbs phenomenon appears near the 
discontinuity; however, the magnitude of the oscillations seems to be independent of 
the mesh size. This is in contrary to conventional FVMs, in which case solutions 
blow up quickly as mesh is refined, unless further stabilization measure (like slope 
limiting or flux limiting) is taken. The last example in ID shows the versatility of 
the proposed FD-FV methods, by solving a scalar nonlinear HCL with non-convex 
flux function. Finally, the first-order operator is extended to solve a 2D Euler 

problem using the framework described in Scction[2j In this 2D example, the superior 
accuracy property proven in Theorem |3.5| is shown to be carried to the 2D case; and 
performance comparison with MUSCL solves show that the proposed second-order 
FD-FV scheme is more efficient than conventional second-order FVMs. 
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RK4. 
RKb 



Mesh size (log scale) Mesh size (log scale) 



Fig. 6.1. Li error of various FD-FV schemes applied to solve \6.1\ : cell-averaged variables 
(left); nodal variables (right) 



For all the tests, the initial condition of cell-averaged value must be computed 
carefully. One choice is to use the exact cell-averaged data obtained from the initial 
condition; an alternative method is to use quadrature points to compute high-order 
accurate cell-averaged values. In particular, if a pth-order FD-FV scheme is tested, the 
initial cell-averaged data is computed to be pth-order accurate approximations to the 
exact value. The former method is adopted for all the tests with discontinuous initial 
conditions; whereas for all tests with smooth initial conditions the second approach 
is employed. 



6.1. One-dimensional advection equations. The first ID advection problem 
has smooth solution and periodic boundary condition: 



ut-f 2u, = 0, {x,t) e [-1, 1] X [0, 1] 

0) = 1 -f ^sin(7ra;), x e [-1, 1] (6.1) 
u{-l,t) = t G [0, 1] 



The exact solution at t = 1 is used to compute the Li errors of both cell-averaged 
variables and nodal variables. Four uniform meshes with number of cells ranging from 
20 to 160 are used to study the convergence behavior. The Li errors are reported in 
logarithmic scales in Figure [O] 
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Fig. 6.2. Li error of various FD-FV schemes applied to solve \6. 2p cell-averaged variables 
(left); nodal variables (right) 



The second problem has smooth solution and Dirichlet boundary condition: 



ut +Mx = 0, {x,t) e [-0.5, 0.5] X [0, 0.5] 



0) 



1 + sin(27ra;) x e [-0.5, 0] 



1 



e [0, 0.5] 



u{-0.5,t) = l+^{t+^fsm{2n{t+^)), i e [0, 0.5] 
u(0.5,^) = 1.0, t e [0, 0.5] 



(6.2) 



Four uniform meshes with number of cells ranging from 20 to 160 are used to study 
the convergence behavior. The Li errors are reported in logarithmic scales in Fig- 
ure lO 



All the simulations use Courant numbers just below the Xmax in Table [4731 I* is clear 
that all the FD-FV schemes are one-order higher than the designed order of , thus 



confirm the superior accuracy property given by Theorem 3.5 



Remark: In order to solve Eqs. (6.2), no special treatment is taken to apply the 
Dirichlet boundary condition at intermediate Runge-Kutta stages. As pointed out 
by Mark H. Carpenter et al. [7], this is not appropriate for high-order Runge-Kutta 
schemes. However, there seems to be no fix of this issue in literature for general 
explicit Runge-Kutta methods. Thus in this test, the fix to RK4 in their work is 
not applied. Nevertheless, the obtained Li convergence behavior is still as expected 
and no degrading of convergence rates is observed. This may be partially due to the 
linearity of the problem. 
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Fig. 6.3. Solution atT = 1.0 obtained by various schemes to solve i 6.3): (1) upper-left, D^'"''- 



RK2; (2) upper-right, Vl'^'^-RKS, 
RK4; (5) bottom-left, D4,«P-i>iased 



(3) middle-left, V, 



3,up — biased 



-RK4; (4) middle-right, T) 



3,up _ 



-RK5; (6) bottom-right, second-order FVM without limiter 



Next, ID advection equation with discontinuous initial condition is solved: 

Ut + Ux = 0, {x, t) e [-1.5, 1.5] X [0, 1] 

r 1 xe [-1.5, -1] U [0, 1.5] 
u{x,Q)=l (6.3) 
[2 X e [-1, 0] 

u(-1.5,t) = u(1.5,t) = 1, te [0, 1] 

Each of the five FD-FV schemes is tested with three uniform meshes with 30, 60 and 
120 cells respectively. These solutions are plotted in Figures [6731 For comparison, so- 
lutions obtained from a second-order FVM using Roe solver, but without any limiting 
strategy, are also presented. 



Since the differential operators "Dx are designed based on Taylor series expansions, 
which are not valid when discontinuity appears, the Gibbs phenomenon is expected 
near the discontinuities. These are clearly observed in Figures |6.3[ However, it is in- 
teresting to see that the magnitude of the spurious oscillations are independent of the 
mesh sizes; whereas in the case of second-order FVM without limiting, the magnitude 
increase slightly when mesh is refined. 
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Fig. 6.4. Li error in density of various FD-FV schemes applied to solve l[6.4y : cell-averaged 
variables (left); nodal variables (right) 



6.2. One-dimensional Euler equations. The first ID nonlinear problem is 
chosen to be the ID Euler equations, in which case the flux function is convex. First, 
ID Euler equations with smooth solutions and periodic boundary condition is tested: 





= 0, 


{x,t) e [-1, 1] X [0, 0.3 


u{x,0) 


= 2 


-f - sin(7ra;j 


p{x,0) 


= 1 


-f - sin(7ra;j 


p{x,0) 


= 1 


-f - sin(7rxj 


w{-l,t) 


— w 


i e [0, 1] 



(6.4) 



Here the conservative vector and flux function are: 



pu 
pu^ + p 
{E+p)ui 



And the system is closed by ideal gas law with specific heat ratio 7 = 1.4: 

P It 



E = 



7-1 



- pu 



To compute the errors, a reference solution is computed by the fifth-order accurate 
FD-FV scheme Pl^.^P-fcmsed.j^j^g ^ uniform mesh with 2560 cells. Each of the five 
schemes is tested using four uniform meshes with number of cells ranging from 20 to 
160. The Li errors in density are reported in logarithmic scales in Figures |6.4[ 



The Li errors in velocity and pressure show similar convergence behavior and are 
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u 


u{w) 


FD-FV scheme 


Cells 


Li errors 


Rates 


Li errors 


Rates 




40 


6.61£;- 3 




4.72£;- 3 






80 
160 


1.80£;-3 
4.70£; - 4 


1.87 
1.94 


1.33£;- 3 
3.50£;-4 


1.83 
1.92 




320 


1.20£;-4 


1.97 


8.98^;- 5 


1.96 




40 


4.05£;-4 




3.54£;-4 




'7^2 WD Tijyo 

Dp '^-tiKd 


80 
160 


5.18£;- 5 
5.88£;- 6 


2.97 
3.14 


4.54£;- 5 
5.11£;- 6 


2.96 
3.15 




320 


8.21£;- 7 


2.84 


7.22E-7 


2.82 




40 


7.66^; - 5 




QME-h 




■7^3 WD— biased Tijy a 


80 
160 


5.20£;-6 
3.36£; - 7 


3.88 
3.95 


AA6E ~ 6 
2.87£;- 7 


3.89 
3.96 




320 


2.13£;-8 


3.98 


1.82£;- 8 


3.98 




40 


4.39£;~ 5 




3.63£;- 5 




I?3."P-RK4 


80 
160 


2.82£;- 6 
1.79£;- 7 


3.96 
3.98 


2.43£;- 6 
1.54£;- 7 


3.90 
3.98 




320 


1.12£;- 8 


3.99 


9.70£; - 8 


3.99 




40 


8.12£;-- 6 




6.63£;- 6 






80 
160 


2.94£; - 7 
9.78£; - 9 


4.79 
4.91 


2.69£;- 7 
9.09E'- 9 


4.62 
4.89 




320 


3.15£:- 10 


4.96 


2.93£:- 10 


4.95 



Table 6.1 



Li error in velocity of various FD-FV schemes applied to solve !j6.4\ 



reported in Table |6.1f]6.2| Note that in this case, the velocity and pressure of cell- 
averaged values are computed by following equations: 

u^u{w)^P^, p^p{w) = {j-1)(e-1^-^) (6.5) 
P \ ^ P J 

The superior accuracy property is confirmed for both cell-averaged values and nodal 
values. 

Next the classical Sod's shock tube problem [3T is tested: the ID Euler equations are 
solved on (a;,i) G [—2, 2] x [0, 0.8] with initial conditions: 

p(a;,0) = 1.0 p(a;,0) = 0.125 

M(a;,0) = 0.0 a;e[-2, 0] , u(x,0)=0.0 a; G [0, 2] 

p(a;,0) = 1.0 p(a;,0)=0.1 

For each of the five FD-FV schemes, three tests are performed using uniform meshes 
with 20, 40 and 80 cells. The solutions of density are plotted in Figures |6.5| Note 
that the solutions are obtained by smaller Courant numbers (reported together with 



the figures) than Xmax in Table 4.3 otherwise the algorithms break down. 



It is again observed that: (1) Gibbs phenomenon appears near the discontinuities; (2) 
the magnitudes of the oscillations around both contact dicontinuity and the shock are 
independent of the mesh sizes. For comparison, solutions obtained using the same 
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P 


p{w) 


FD-FV scheme 


Cells 


Li errors 


Rates 


Li errors 


Rates 




40 


7.54£; - 3 




5.98^;- 3 






80 
160 


2.01£;-3 
5.17£;-4 


1.91 
1.96 


1.61£;- 3 
4.14£;-4 


1.90 
1.96 




320 


1.31£;-4 


1.98 


1.05£;-4 


1.98 




40 


4.62£;-4 




4.10£;-4 






80 
160 


6.68£;- 6 


2.96 
3.15 


5.23E'- 5 
5.83£;- 6 


2.97 
3.17 




320 


9.33£;- 7 


2.84 


8.26£;- 7 


2.82 




40 


7.85^;- 5 




6.49£;- 5 




■j^S, up —biased p^J^^ 


80 
160 


5.29£;-6 
3.39£; - 7 


3.89 
3.96 


4.26£;- 6 
2.72£;- 7 


3.93 
3.97 




320 


2.13£;-8 


3.99 


1.71£;- 8 


3.99 




40 


4.32£;- 5 




3.27£;- 5 






80 
160 


2.82£;- 6 
1.78£;- 7 


3.94 
3.98 


2.21£;- 6 
1.41£;- 7 


3.88 
3.97 




320 


1.12£;- 8 


3.99 


8.89£;- 8 


3.99 




40 


8.40£; - 6 




6.51E'- 6 






80 
160 


3.09£; - 7 
1.02£;-8 


4.76 
4.92 


2.64£; - 7 
9.13£;- 9 


4.62 
4.86 




320 


3.27£:- 10 


4.97 


2.94£:- 10 


4.96 



Table 6.2 



Li error in pressure of various FD-FV schemes applied to solve \6.Ji\ 



meshes by a second-order FVM using Roe flux and without any limiting strategy are 



also provided in Figure 6.5 Note that the unlimited FVM performs poorly near the 
shock: no data is obtained at all on the finest mesh. Hence it is observed, and conjec- 
tured that the proposed FD-FV methods have some inherent stabilization mechanism 
that is not possessed by conventional FVM methods. 

It must be emphasized, however, the author do not mean by this observation, that the 
present FD-FV schemes are competitive to modern FVMs, especially when disconti- 
nuity is expected. The reason is that nonlinear stability analysis of FVMs are very 
mature, and many techniques are developed over decades to ensure certain stability 
property (like TVD, TVB, or discrete maximum principle) for finite volume schemes. 
Nevertheless, it worths to study the nonlinear stability property of the proposed FD- 
FV schemes to enhance the robustness of the methods; this is an important direction 
of future work on FD-FV schemes. 

The last test with ID Euler equations is a more challenging one, the Shu-Osher 
problem [33] . This problem solves flow obtained by a moving Mach 3 shock wave in- 
teracting with a sinusoidal density proflle. The computational domain is a; G [—5, 5], 
with initial conditions: 

p{x, 0) = 3.857143 p{x, 0) = 1 + i sin(5a;) 

M(a;,0) = 2.629369 x e [-5, -4] , u\x,Q) ^ QX) ' x G [-4, 5] 

p{x, 0) = 10.33333 p{x, 0) = 1.0 

The problem is solved until T = 1.8. In this case, not all the five FD-FV schemes 
succeed in producing solutions. In particular, density profiles obtained by I?^^"''-RK2, 
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Fig. 6.5. Density at T = 0.8 obtained by various schemes applied to solve Sod's problem: 
(1) upper-left, Vl'"^^ -RK2, CFL = 0.8; (2) upper-right, V'^'"^-RK3, CFL = 0.2; (3) middle- 
left, Vl'''''~^^°'"^'^-RK4, CFL = 0.4; (4) middle-right, Vl'"^-RK4, CFL = 0.1; (5) bottom-left, 
^4,up-biased_j^j^^^ CFL = 0.2; (6) bottom-right, second-order FVM without limiter, CFL = 0.8 




Fig. 6.6. Density profiles of Shu-Osher problem obtained by MUSCL-RK2 (500 cells) and 
FD-FV schemes (250 cells) at T = 1.8; global view (left); local view (right) 



X>2^"P-RK3 and X>3,«p- 



-biased 



RK4 with a reduced Courant number are presented in 



Figures 6.6 6.7 For comparison, the density profile obtained by the MUSCL scheme 
with RK2 and van Albada hmiter |38j is also plotted. In the same figure, the number 
of cells for FVM is twice as large as number of cells for FD-FV schemes, so that the 
numbers of DOF are the same. The reference solution is computed using MUSCL- 
RK2 with 10000 uniform cells. 
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Fig. 6.7. Density profiles of Shu-Osher problem obtained by MUSCL-RK2 (1000 cells) and 
FD-FV schemes (500 cells) at T = 1.8; global view (left); local view (right) 



Comparing the two second-order schemes (MUSCL-RK2 and 2?;|.'"P-RK2), given the 
same number of DOF: (1) the FD-FV scheme provides much more accurate magnitude 
of the solution, indicating the proposed scheme is much less dissipative than FVM 
with the same order of accuracy; (2) on the other hand, the phase error obtained 
from the second-order FD-FV scheme is larger than MUSCL-RK2. This phase error 
is greatly reduced if higher-order FD-FV methods are used. 

6.3. One-dimensional scalar conservation law with non-convex flux. 

The last ID test case concerns non-convex flux [17j . in which case approximates 
Riemann solver in FVMs is much harder to design. The governing equation is given 
by the scalar HCL: 

ut + fiu),=0, f{u) = ^{u^ -l)iu^ -4) (6.6) 

The problem is solved on {x,t) £ [—2, 2] x [0, 0.04] with initial condition: 

it(a;,0) = -3, xe[~2,0]; u{x,0)^3, x e [0, 2] 
The exact solution at T = 0.04 is given by: 



u{x, T) 




X < -19.5r 
-19. 5r <x <Q 
< X < 19. 5r 
19. 5r < X 



Here g{ ) satisfies u = f'{g{u)). All the five FD-FV schemes are tested. The first 
series of results are obtained using 80 uniform cells, and the results are plotted with 
the exact solution in Figures 6 



Note that the magnitudes of the oscillations are not symmetric about x — Q, this is 
due to the fact that x = Q coincides with a grid point, and the initial condition at 
this point is set to be equal to 3, i.e. the right status. If odd number of uniform cells 



is used, for example 81 cells, the asymmetry disappears, as shown in Figures 6.9 



It is also interesting to observe that, using odd number of cells, the spurious oscillations 
near the weak discontinuities almost disappear. 

6.4. Two-dimensional Euler equations. Finally, a 2D problem, the isen- 
tropic vortex advection problem |32) . is tested to illustrate the superior accuracy 
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Fig. 6.8. Solution to the scalar conservation law with non-convex flux with 80 cells using: 
Vl''''^-RK2, Vl'^'^-RKS and V't'^P''''"""' -RK4 (left); Vl'""^ -RK4 and ©^'"''"""'"^''-RiCS (right) 





y' 

J 







Fig. 6.9. Solution to the scalar conservation law with non-convex flux with 81 cells using: 

Vi'''P-RK2, Vl'^'P-RKS and JJ^^^P-bio.sed ^^^j^^. J)3,up ^4,up-i,ia.ed_^^^ ^^.^j^^j 



property of the proposed FD-FV scheme. In particular, I?^'"^-RK2 is employed, and 
its performance is compared with conventional second-order FVM, i.e. MUSCL~RK2, 



using van Albada slope limiter. Note that the operators and Dy in Eqs. (2.13 



2.14) are chosen to be first-order accurate upwind finite difference operators. The 



governing equation is the 2D Euler equations: 



wt + f{w)x + giw)y = 



pu 
pv 



( 



pu 
pu^ + p 



\ 



puv 
\{E + p)uJ 



g{w) 



pv 
puv 
pv^ + p 
\{E+p)vJ 



(6.7) 



This system is closed by: 



E=^ + Ip{u^ + v^), 7 = 1.4 
7 — i z 



The computation is performed on {x,y,t) = [—5, 5] x [—5, 5] x [0, 10], with periodic 
boundary condition on all four edges. The initial condition is an isentropic vortex 
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Grid sizes 


Li errors 


Rates 


Li errors 


Rates 


P 


P 


20 X 20 


1.94£; + 




1.87£; + 




40 X 40 


8.48£; - 1 


1.19 


8.12£;- 1 


1.20 


80 X 80 


2.56£;- 1 


1.73 


2.41£;- 1 


1.75 


160 X 160 


6.98£;- 2 


1.88 


6.26£; - 2 


1.95 




u 


u{w) 


20 X 20 


6.12£; + 




6.36E + 




40 X 40 


2.74£; + 


1.16 


2.79£; + 


1.19 


80 X 80 


8.25£;- 1 


1.73 


8.19^;- 1 


1.77 


160 X 160 


2.18E;- 1 


1.92 


2.14£;- 1 


1.94 




V 


v{w) 


20 X 20 


5.46£; + 




5.76E + 




40 X 40 


2.06£; + 


1.40 


2mE + 


1.47 


80 X 80 


5.92E- 1 


1.80 


5.83E - 1 


1.84 


160 X 160 


1.55E- 1 


1.93 


1.50£;- 1 


1.97 




P 


p{w) 


20 X 20 


2A6E + 




2.41£; + 




40 X 40 


1.09£; + 


1.17 


1.03£; + 


1.23 


80 X 80 


3.26£;- 1 


1.74 


2.91E- 1 


1.82 


160 X 160 


8.38E; - 2 


1.89 


7.56£; - 2 


1.95 



Table 6.3 

Li error in w and w of T)^'^'' -RK2 applied to solve isentropic vortex advection problem 



given by: 



w(a;,2/,0) = 1- — exp |^-(l-a; -y 
v{x,y,0) = 1 + ;^exp ( ;^(1 - - y'^] 

ZTT y / 

p{x,y,0) = (l - ^l'^ (l~x^-y 



2^ 



p{x, y, 0) = (^1 - ^^^^11'^ cxp (1 - x2 - y^] 

This initial condition leads to a vortex with uniform entropy p/ — 1.0. The constant 
e is chosen to be e = 5, and the exact solution at one period T = 10 equals the initial 
condition. Four uniform meshes are employed to test the convergence behavior with 
grid sizes ranging from 20 x 20 to 160 x 160. The computed Li errors are reported in 
Table [Q 



It is clear from the table that the asymptotic convergence rates in all components 
are second-order. Note that the coarsest mesh corresponds to using only 8 cells per 
dimension to resolve the vortex; hence the computed Li errors are large. The com- 
putational costs of this FD-FV scheme is compared to MUSCL-RK2, which has the 
same order of accuracy in both space and time. To this end, Table |6.4| presents the 
comparison results obtained by choosing the same mesh size for both schemes, as 
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Scheme 


Mesh size 


Ndof 


Ep 


Niter 


^cpu (SGC.) 


FD-FV 


20 X 20 


4,960 


1.87i;-h0 


200 


20 


FVM 


20 X 20 


1,600 


2.23E + 


200 


18 


FVM 


35 X 35 


4,900 


1.26E + 


350 


92 


FD-FV 


40 X 40 


19,520 


8.12E~ 1 


400 


158 


FVM 


40 X 40 


6,400 


1.02£; + 


400 


139 


FVM 


70 X 70 


19, 600 


3.76E - 1 


700 


731 


FD-FV 


80 X 80 


77, 440 


2.41£;- 1 


800 


1,252 


FVM 


80 X 80 


25,600 


2.97^; - 1 


800 


1,101 


FVM 


140 X 140 


78,400 


1.25E- 1 


1,400 


5,852 


FD-FV 


160 X 160 


308, 380 


6.26£; - 2 


1,600 


9,945 


FVM 


160 X 160 


102,400 


1.07^;- 1 


1,600 


8,730 


FVM 


278 X 278 


309, 136 


4.73E - 2 


2,780 


45,938 



Table 6.4 

Performance comparison between second-order FD-FV fDx'^^ -RK2) and second-order FVM 
(MUSCL-RK2 with van Albada lirniter): N^^f ^ number of DOFs; Ep - Li error in density for 
cell-averaged data; Niter - number of iterations; Tcpu - computational time in seconds 



well as the comparison results using roughly the same number of DOFs. Here all the 
tests are performed using C-|— |- programs and g-|— I- compiler under 64bit CentOS 5.8 
system with single processor (Intel(R) Xeon (TM) CPU 3.00GHz). 

It can be concluded from the table that: (1) given the same mesh size, although the 
number of DOFs of pi'"P-RK2 is two times larger than MUSCL-RK2, the CPU times 
arc almost the same; and the FD-FV scheme gives much better results than FVM; 
(2) given the same number of DOFs, X'^'"^'-RK2 provides less accurate solutions than 
MUSCL-RK2, but the former approaches asymptotic convergence faster than the lat- 
ter as the meshes arc refined - the difference between the Li errors is less significant 
when meshes are refined; (3) given the same number of DOFs, I?^'"*'-RK2 is much 
faster than MUSCL-RK2. 

Remark: This comparative study only considers explicit time-integrators, in which 
case the cost of FVMs is mainly due to numerical flux evaluations. However, implicit 
methods are often employed in practice to allow larger time steps. When implicit 
time-integrators are used, a major additional computational cost is from the nonlin- 
ear solver, for example, by Newton methods. In this case, the connectivity plays an 
important role in determining the efficiency of the algorithm. Thus the statements 
made here may or may not applies to implicit methods, which will be investigated in 
future work. 

7. Concluding remarks. A hybrid finite difference-finite volume approach is 
proposed for solving general first-order hyperbolic conservation laws. It is a departure 
from conventional finite difference schemes and finite volume schemes in that both cell- 
averaged values and nodal values are counted as dependent variables and evolved in 
time. By distinguishing between cell-averaged DOFs and nodal DOFs, the proposed 
FD-FV schemes are: (1) natiiral conservative like FVMs, but work for general flux 
functions and do not require any exact or approximated Ricmann solver; (2) easy 
to extend to high-order accuracy in space like FDMs. The present work focuses 
on discretization in space, and integration in time is achieved by method of lines; in 
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particular, explicit Runge-Kutta methods are adopted for simplicity. It is theoretically 
proven and numerically illustrated that the proposed FD-FV methods have superior 
spatial accuracy property than conventional FDMs and FVMs. This means that the 
achieved order of accuracy in space is one-order higher than the designed order of 
accuracy of the discrete differential operator applied in the FD-FV scheme. 
When shock occurs in the solutions, second-order FVM fail if no other stabilization 
mechanism is applied, especially when the mesh is sufficiently refined. However, it 
is observed for the proposed FD-FV schemes, that the magnitudes of the oscillations 
near both contact dicontinuties and shocks are independent of the mesh size. Hence it 
is conjectured that the FD-FV schemes have some inherently stabilization mechanism 
that worths further investigation. 

There are remaining issues, however, to be addressed in future work. For example, 
nonlinear stability of the methods is to be investigated, and the present schemes need 
to be enhanced to avoid any spurious oscillations near discontinuities. In addition, 
since the conservation law alone does not prevent non-physical solutions, entropy 
inequality should be considered in designing practical FD-FV schemes. 
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